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ABSTRACT 


A  theoretical  method  for  representing  the  physical  characteristics 
of  a  portion  of  a  mechanical  piping  system  in  terms  of  a  stiffness 
matrix  involving  only  those  joints  at  which  the  sub-system  is  connected 
into  the  overall  system  is  presented  herein.   The  use  of  this  method  in 
reducing  the  size  of  the  matrices  that  are  needed  for  finding  the  normal 
vibration  frequencies  of  a  piping  system  is  discussed,  and  an  outline  for 
a  digital  computer  program  to  use  this  method  is  presented. 
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CHAPTER  I 
INTRODUCTION 

1.1  General  Remarks 

There  are  many  aspects  to  be  considered  in  the  design  of  any  pip- 
ing system;  stresses,  strains,  pressure  drops,  and  corrosion,  to  name 
a  few.  This  thesis  deals  with  only  one  aspect  of  the  overall  problem, 
that  of  the  dynamic  behaviour  of  the  system  as  represented  by  its  nat- 
ural frequencies  of  vibration,  and  presents  a  method  whereby  a  large 
system  may  be  systematically  attacked  by  considering  small  sub-systems 
in  order  to  improve  computer  utilization. 

Except  for  very  simple  or  trivial  problems  the  number  of  calcula- 
tions involved  and  the  size  of  the  matrices  used  make  the  use  of  a 
large  digital  computer  mandatory.   No  digital  computer  program  has  been 
written,  however,  but  an  outlined  plan  of  approach  is  included  in  this 
thesis. 

1.2  Scope  of  Work  Presented 

Baird  [1]*  has  presented  a  method  for  determining  analytically  the 
undamped  natural  frequencies  and  mode  shapes  of  large  three  dimensional 
piping  systems  by  use  of  stiffness  matrices.   He  also  has  made  use  of 
topological  matrices  to  express  mathematically  certain  physical  or 
boundary  conditions.  The  most  obvious  inherent  disadvantage  of  this 
method,  however,  is  the  large  size  of  the  matrices  which  may  be  re- 
quired with  the  result  that  supplementary  memory  of  some  sort  must  be 


♦Numbers  appearing  in  this  manner  refer  to  the  bibliography,  page  41, 


used  in  the  computer  solutions,  thus  increasing  the  computation  time 
considerably. 

A  technique  is  presented  in  this  thesis  whereby  small  portions  of 
the  system,  called  ganglia  or  sub-systems,  can  be  analyzed  and  reduced 
to  some  sort  of  pseudo-elements  having  certain  calculated  characteris- 
tics with  the  result  that  only  the  reactions  at  the  ends  of  these  pseudo- 
elements  are  of  interest,  and  all  the  internal  characteristics  are  en- 
gulfed in  the  general  stiffness  matrix  for  each  ganglion.   These  individ- 
ual sub-systems  can  be  used  to  build  up  an  entire  system  whose  vibra- 
tional characteristics  can  be  represented  far  more  compactly  than  can 
those  of  the  unsimplified  system. 

An  outline  for  a  proposed  computer  program  to  utilize  this  method 
is  included  in  Chapter  IV. 

1.3  Definitions 

A  number  of  definitions  must  be  established  and  the  reader  is  warn- 
ed that  they  are  somewhat  at  variance  with  those  of  Baird  [1],  The 
nomenclature  used  (see  section  1.4)  will  also  be  at  variance  with  Baird' s. 
Element  -  A  piece  of  the  system  having  two,  and  only  two,  ends  and 
no  internal  branches.   Normally  an  element  will  be  a  simple  length 
of  pipe.  A  tee,  though,  which  is  a  more  complicated  item,  should 
be  thought  of  as  three  elements  joined  at  a  common  point.   Other 
items  of  piping  hardware  may  be  analyzed  similarly. 
Universe  -  This  term  is  used  in  the  same  sense  in  which  it  is  used 
in  statistics,  namely  the  larger  body  of  effects  from  which  the 
system  is  taken.   In  this  case  the  universe  for  a  piping  system 
consists  of  all  the  boundary  forces  and  constraints,  foundations, 
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etc.  which  affect  the  system.  When  a  sub-system  is  isolated,  its 
universe  is  the  larger  system  itself. 

Node  -  A  terminal  point  of  an  element  is  called  a  node.   It  is 
noted  that  a  node  is  a  theoretical  point  and  not  a  physical  item. 
There  must  be  at  least  one  element  end  at  each  node,  though  there 
may  be  any  number  greater  than  one.  Nodes  may  be  further  classi- 
fied as: 

1.  Trivial  nodes  -  Those  nodes  upon  which  only  two  element 
ends  are  incident. 

2.  Boundary  Nodes  -  Those  nodes  at  which  there  are  forces  and 
constraints  applied  to  the  system  by  the  universe.  For  pur- 
poses of  convenience  when  dealing  with  coordinate  transforma- 
tions it  is  required  that  only  one  element  be  incident  upon  a 
boundary  node.   In  many  cases  this  may  require  the  introduc- 
tion of  an  element  connecting  a  primary  node  to  the  boundary 
node  at  which  the  system  connects  to  the  universe.  These 
primary  and  boundary  nodes  may  have  the  same  geometrical  loca- 
tion, though  they  are  topologically  distinct. 

3.  Primary  Nodes  -  All  nodes  of  the  system  that  are  not  trivial 
or  boundary  nodes. 

Path  -  A  sequence  of  elements  connecting  two  nodes  such  that  there 

are  only  trivial  nodes  connecting  these  elements. 

Primary  Path  -  A  path  connecting  two  primary  nodes  or  a  primary  node 

and  a  boundary  node.   A  primary  path  may  contain  only  trivial  nodes 

on  its  interior. 

Primary  Element  -  An  element  made  up  of  all  the  elements  on  a  primary 

path. 


State  Vector  -  A  column  matrix  of  displacements  and  corresponding 
forces  at  any  point.   Since  in  general  it  requires  three  translations 
and  three  rotations  to  describe  motion  completely  the  state  vector 
is  a  12  x  1  matrix.   For  the  purposes  of  this  thesis  it  has  been 
convenient  to  make  the  first  six  elements  of  the  column  the  displace- 
ments in  the  six  independent  degrees  of  freedom,  and  the  last  six 
elements  the  corresponding  forces. 

Transfer  Matrix  -  A  square  (12  x  12)  matrix  representing  the  physi- 
cal characteristics  of  an  element  such  that  the  state  vector  at  one 
end  will  equal  the  product  of  the  transfer  matrix  and  the  state 
vector  at  the  other  end,  in  that  order. 

Dynamic  Transfer  Matrix  -  A  transfer  matrix  that  takes  into  account 
the  mass  of  the  elements  and  the  frequency  of  vibration  in  addi- 
tion to  the  static  properties  (i.e.,  Young's  modulus,  Poisson's 
ratio,  second  moments  of  area,  etc.). 

Ganglion  -  A  sub -system  "removed"  from  the  system  at  primary  nodes. 
"Removed"  is  used  to  indicate  that  the  ganglion  is  considered  as  a 
separate  system  whose  boundary  nodes  are  the  nodes  at  which  it  was 
connected  into  the  overall  system. 

1.4  Nomenclature 

A  capital  letter  will  normally  be  used  to  represent  a  matrix  and 

appropriate  sub-  and  super-scripts  may  be  used  to  further  define  it  as 

follows: 

a  A  C 
b  A  d 

a  -  number  of  rows 

b  -  number  of  columns 


c  -  may  be  1)  blank,  2)  T  indicating  the  transpose  of  the 
matrix,  or  3)  -1  indicating  the  inverse  of  the  matrix 
d  -  used  for  identification  and  described  as  required  in 
context. 
Brackets,  [  ],  may  be  used  to  delineate  a  matrix  or  the  partition  of  a 
matrix  if  the  meaning  is  thereby  clarified.  They  will  always  be  used 
when  the  elements  of  a  matrix  are  exhibited  in  extensio.  Braces. 


,[}. 


will  be  used  to  permit  the  writing  of  the  terms  of  a  column  vector  as 
a  row  vector  to  save  space.   The  term  Diag <    >  will  mean  that  the  n 
terms  inside  the  braces  are  the  n  diagonal  terms  of  an  otherwise  null 
matrix,  though  the  terms  themselves  may  represent  square  sub-matrices. 

In  no  case  will  more  than  one  capital  letter  be  used  for  a  matrix; 
thus  the  expression  AB  will  mean  the  product  of  matrices  A  and  B  in  that 
order  and  will  imply  that  A  is  conformable  to  B  for  multiplication.   One 
should  recall  that  matrix  multiplication  is  distributive  and  associative 
though  not  necessarily  commutative. 

All  nodes  will  be  numbered,  and  for  convenience  the  boundary  nodes 
will  be  given  the  lowest  numbers.   The  ends  of  the  elements  will  be 
identified  by  lower  case  English  letters,  two  per  element. 

Notation  generally  follows  that  of  Pestel  and  Leckie  [2].  The 
following  are  the  matrix  symbols  used. 

A    Node  incidence  matrix 

D,  F  Displacement  and  force  vectors  (or  partitions)  respectively 

G    Coordinate  rotation  matrix 

K    Dynamic  stiffness  matrix 

K^   Node  stiffness  Matrix 
Jicj   Frequency  determinant 


K    Primitive  stiffness  matrix 
P 

T    Transformation  matrix 
U    Transfer  matrix 
Z    State  vector 
By  the  use  of  subscripts  any  of  the  above  symbols  may  be  used  to 
represent  an  element  or  partition  of  the  full  matrix,  thus  Z.  -  would 
be  an  element  or  partition  (depending  upon  the  context)  of  matrix  Z, 
namely  the  one  appearing  in  row  1,  column  2. 
The  following  scalar  quantities  are  used: 

B    The  number  of  boundary  nodes  of  a  system 
M/2  The  number  of  primary  paths  in  a  system 

(Therefore  M  is  the  number  of  primary  element  ends  in  the 
system.) 
N    The  number  of  nodes  (boundary  and  primary)  in  a  system 
The  use  of  a  prime,  ('),  will  generally  indicate  that  the  quantity  desig- 
nated refers  to  a  sub-system. 
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CHAPTER  II 
THE  TRANSFER  AND  STIFFNESS  MATRIX  SOLUTION 

A  brief  summary  of  the  solution  by  use  of  stiffness  matrices  for 
the  characteristic  (normal)  frequencies  of  a  piping  system  with  general 
topology  is  given  in  this  chapter  in  order  to  refresh  the  reader's  memory. 
For  further  details  of  the  theory  see  Baird  [1],  but  again,  be  aware  of 
notational  differences. 

Before  we  discuss  the  dynamic  stiffness  matrix  approach  as  develop- 
ed in  fl]  and  [2]  it  is  necessary  to  remind  the  reader  that  this  direct 
mathematical  approach  deviates  somewhat  from  the  approach  to  planning  a 
digital  computer  program  because,  although  on  paper  we  can  manipulate 
symbols,  in  the  computer  we  must  deal  with  numbers.   It  is  assumed  that 
the  user  of  this  method  has  available  the  necessary  dynamic  transfer 
matrices  for  whatever  fundamental  elements  are  involved  and  can  compute 
(in  one  way  or  another)  the  geometrical  matrices  necessary  to  relate  all 
quantities  to  some  single  coordinate  system. 

a.  Calculate  a  transfer  matrix  for  each  primary  element,  eliminat- 
ing the  trivial  nodes  along  the  path.   Thus: 

U.    -  U   G  U  .   G  -  ...U,  0G-U,  , 
k,m    n,m  n  n-l,n  n-1     1,2  1  k,l 

where  k  and  m  are  element  ends  at  primary  and/or  boundary  nodes,  and  n, 

n-1,  ...1  are  the  trivial  nodes  along  the  path.   G  ,  G  , ,...G-  are  the 

n   n-1     1 

coordinate  rotation  matrices  providing  the  necessary  coordinate  align- 
ment at  each  node.   The  resultant  transfer  matrix,  U,    is  in  the  coor- 

k,m 

dinate  system  of  the  element  end  at  the  node  at  which  the  path  starts, 
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end  k  in  this  case.   The  relationship  between  k  and  m  is  now: 


Z  -  U.    Z, 

m    k,m  k 


2.1 


b.  Transform  the  transfer  matrix  of  each  primary  element  into  a 
dynamic  stiffness  matrix  as  indicated  in  what  follows.  The  user  is  cau- 
tioned to  establish  and  adhere  to  an  order  within  his  state  vector. 
The  order  used  by  Baird  [1]  is  recommended.   Referenced  to  his  local 
coordinate  system  it  is: 


Z  - 


L 


Displacement  in  an  x  direction 
Displacement  in  a  y  direction 
Displacement  in  a  z  direction 
Rotation  about  an  x  axis 
Rotation  about  a  y  axis 
Rotation  about  a  z  axis 
Force  in  an  x  direction 
Force  in  a  y  direction 
Force  in  a  z  direction 
Moment  about  an  x  axis 
Moment  about  a  y  axis 
Moment  about  a  z  axis 


Partitioning  the  state  vectors  into  force  and  displacement  vectors, 
and  the  dynamic  transfer  matrix  into  four  6x6  matrices  yields: 


m 


m 


U 


11 


u 


21 


U 


12 


U 


22 


D, 


2.2 


An  explanation  of  the  use  of  the  lower  case  f  will  follow  shortly.  The 
above  equations  then  can  be  written: 


12 


in 


unDk  +  ui2fk 


m 


U21Dk   +   U22fk 


Solving  the  first  of  these  for  f,  gives 

f,    -    U.'J  D   -  U.'i  U-.  D. 
k         12  m      12   11  k 


and  substituting  this  into  the  second  gives 
F 


m 


D21Dk     +     °22<Di2Dm     "     UU°llV 


Reassembling  these  last  two  equations  in  matrix  form,  we  have 


-fk 

■ 

"m  un 

-ll      ' 

F 
m 

U21-U22Di2Ull 

U22°U 

.  



m 


2.3 


There  is  a  sign  convention  that  must  be  carefully  delineated  at  this 
time.   In  equation  2.1  the  force  elements  of  the  state  vector,  Z  ,  (shown 
in  2.2  as  the  partitioned  sub-matrix  f  )  are  the  forces  applied  to 
element  km  by  the  node  at  k,  whereas  the  force  elements  of  the  state 
vector  Z  are  the  forces  applied  by_  element  km  to  the  node  at  m.   It 
is  obvious  then  that  the  negatives  of  the  elements  of  sub-matrix  f, 
represent  the  forces  applied  by_  element  km  to  the  node  at  k;  and  if  we 
define  a  new  vector  F,  ■  -f,  ,  then  both  F,  and  F  represent  forces  ap- 

K        HI  K         m 

plied  by  the  element  to  the  nodes.  Then  equation  2.3  can  be  rewritten: 


i-  mJ 


V 


m 


D, 


D 

!—  m. 


2.4 


where  K,   ,  the  dynamic  stiffness  matrix  for  the  element  km,  is  the  square 
Tc,m 

matrix  in  2.3. 
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c.  The  dynamic  stiffness  matrix  of  each  primary  element  is  now 
known  in  the  coordinate  system  of  the  element  end  at  which  calculations 
started.   In  order  to  relate  all  these  matrices  to  each  other  it  is 
necessary  that  some  global  coordinate  system  and  reference  axes  be 
chosen  and  that  all  the  dynamic  stiffness  matrices  be  transformed  into 
this  system.   Detailed  discussion  of  the  necessary  geometric  matrices 
and  multiplications  may  be  found  in  [1]  and  [3]  and  will  not  be  repeated 
here. 

d.  The  system  now  consists  of: 

M/2  primary  nodes 

B   boundary  nodes 

N   nodes  (boundary  and  primary) 

zero  trivial  nodes 

e.  A  topological  matrix,  A,  relating  the  incidence  of  element  ends 
on  the  nodes,  is  now  generated  as  follows.   Each  of  its  elements  is  in 
actuality  a  6  x  6  matrix,  0  or  I,  0  representing  a  6  x  6  null  matrix  and 

I  representing  a  6  x  6  identity  matrix.  The  topological  matrix  will  be 

T 
generated  in  its  transposed  form,  namely  A  ,  since  both  forms  are  needed 

T 
in  subsequent  calculations.  A  will  be  6N  x  6M  in  size.  The  following 

rules  obtain: 

T 
A        will  be  I  if  end  j  is  incident  on  node  i. 
i » J 

T 

A  .  .     will  be  0  if  end  i  is  not  incident  on  node  i. 
i,J 

I  and  0  are  6x6  matrices  as  defined  above.   To  permit  doing  this,  a 
certain  order  must  be  established  both  for  the  nodes  and  the  element 
ends,  and  this  order  will  have  to  be  adhered  to  for  subsequent  assemblages 
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and  calculations.   It  will  be  helpful  if  the  low  numbered  nodes  (i.e., 

T 
the  upper  rows  of  the  matrix  A  )  are  the  boundary  nodes  as  this  will 

T 
minimize  subsequent  reordering.  As  an  example,  A  is  shown  for  the 

simple  system  given  below. 


i  la 


h  \ 


Nodes  are  numbered, element  ends  are  lettered.  Nodes  1  and  2  are  boundary 
nodes,  nodes  3  and  4  are  primary  nodes. 


=  1 
2 
3 


£ 


jibcdefgh^ 
10000000 

00000001 

01101000 

OOOIOIIO 


Row  of  all  element  ends 


Column  of  nodes,   note   the  boundary  nodes   listed   first 

T 

A  is  the  quantity  within  the  brackets.   The  row  of  element  ends  and 

the  column  of  nodes  are  shown  only  to  indicate  how  the  matrix  is  generat- 
ed. 

T   T 

As  implied  herein,  A  is  [A  ] 

f.  All  the  primary  element  dynamic  stiffness  matrices  are  now 

assembled  into  a  large  matrix  called  the  primitive  stiffness  matrix,  K  . 

P 

T 
It  is  necessary  that  the  order  of  assembly  correspond  to  that  of  A  . 
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2.5 


Defining  <  F  F,  . . . >  as  F  (subscript  E  referring  to  ends  of  elements) 

and  I  D  D,  . . .  I 
La  b    J 


as  D   then 


E      p  E 

where  K  is  the  large  square  matrix  in  equation  2.5. 
P 


2.6 


g.  The  sura  of  the  forces  exerted  on  any  node  by  the  element  ends  in- 
cident upon  it  must  be  zero  (except  for  boundary  nodes  where  the  constrain- 
ing forces  are  the  "closing"  forces)  and  the  displacements  of  all  the  ends 

incident  upon  any  given  node  must  be  identical.  The  topological  matrices 

T 
A  and  A  permit  these  conditions  to  be  expressed  mathematically.  Using 

the  subscript  N  to  refer  to  nodes*  just  as  subscript  E  refers  to  ends 

of  elements, 


and 


FN  "  A  FE 


DE  "  AV 


F„  is  a  column  vector  of  the  forces  exerted  on  the  nodes  and  D„  is  a 

N  N 

column  vector  of  the  node  displacements.  Therefore: 


F   =  A  F 
N        E 

=  AT  (K  D  ) 
P  E 

-  AT  K  A  DM. 
p    N 
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T 
Defining  A  k  A  as  the  node  stiffness  matrix,  K^,  we  then  have  the  re- 
lationship between  node  forces  and  displacements: 


N 


h\ 


2.7 


h.   P„  is  a  column  of  N  6x1  matrices  and  all  but  B  of  these  will 
N 

have  all  the  elements  identically  equal  to  zero.  F  ,  K^  and  D  are  now 

reordered  so  that  the  B  uppermost  6x1  matrices  of  F„  and  D„  are  the  force 

N      N 

and  displacement  vectors  respectively  of  the  B  boundary  nodes.  This  step 
is  unnecessary  if  the  proper  order  has  been  previously  established  as  in- 
dicated in  2.e.  Partitioning  these  matrices  now  yields: 


6Bp 
1  rNB 


6(N-B), 


^11 
*N21 


12 


22 


NB 


N2 


■;•- 


i.   Since  boundary  constraints  are  known  in  the  local  coordinate 
systems  of  the  various  boundary  nodes  we  now  transform  back  to  local 
coordinates.  Up  to  this  point  all  reordering  has  been  done  using  entire 
6x1  or  6x6  submatrices.  This  has  been  done  in  order  to  preserve  the  in- 
tegrity of  the  partitions  of  the  state  vector  so  that  we  can  make  the 
transformation  back  to  local  coordinates  by  using  the  same  transforma- 
tion matrices  generated  earlier  to  transform  the  system  to  global  coordi- 
nates. Now  with  F   and  D   transformed  to  local  coordinates  we  may  find, 

NB       NB 

upon  examination  of  the  boundary  conditions,  that  some  of  the  elements  of 
F   are  also  equal  to  zero.  Reordering  the  matrices  again,  this  time 
element  by  element  instead  of  by  submatrices,  so  that  the  non-zero  ele- 
ments of  F%m  and  the  corresponding  zero  elements  of  Dim  occupy  the  lead- 
NB  NB 

ing  positions  in  their  respective  vectors,  and  then  partitioning  the 
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matrices  we  find: 


NB 


*  * 

Sll  ^12 

*  * 
*N21  ^22 


N2 


2.8 


NB 


is  a  Pxl  vector  where  P  is  the  number  of  non-zero  boundary  condi- 


tions and  is  less  than  or  equal  to  6B.  The  asterisk  is  used  to  in- 
dicate local  coordinates  and  the  reordering  indicated  in  the  text. 
From  equation  2.8  it  can  be  seen  that: 


* 

* 

* 

0 

= 

* 

0 

+ 

*N22 

DN2 

Since  the  elements  of  D.TO  are  not  zero,  and  are  independent  it  appears 
*  N2 


that  Kw^  must  De  singular,  i.e., 


Defining 


*N22 
* 

*N22 


2.9 


as  the  frequency  determinant,  and  calling  it 


■S, 


for  convenience,  we  can  see  that  it  will  be  of  order  (6N-P)  x  (6N-P). 


Since  some  of  the  elements  of   K^   are  frequency  dependent  it  can  be 
seen  that  there  will  be  some  values  of  frequency  for  which  the  frequency 
determinant  becomes  zero.  These  frequencies  will  be  the  natural  fre- 
quencies of  the  system.  We  are  unable  to  solve  for  these  frequencies 
directly  since  it  would  not  be  feasible  to  obtain  a  literal  expression, 
ultimately  expressible  as  a  polynomial  (or  worse)  in  the  unknown  fre- 
quencies.  Therefore  we  calculate  the  value  of  the  frequency  determin- 
ant for  many  different  values  of  frequency  (using  a  high  speed  digital 
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computer)  and  plot  the  value  of  the  determinant  versus  the  frequency. 
From  this  plot  we  should  be  able  to  pick,  at  least  approximately,  the 
natural  frequencies.   It  is  then  possible  to  refine  these  answers  by 
making  additional  calculations  over  much  smaller  intervals  and  plotting 
these  results  on  an  expanded  scale.   If  a  lumped  mass  model  of  the  system 
has  been  used  throughout,  there  will  be  a  finite  number  of  natural  fre- 
quencies.  If  however,  a  distributed  mass  model  has  been  used  then  there 
will  be  an  infinite  number  of  natural  frequencies. 

j.   It  can  easily  be  seen  that  the  sizes  of  the  matrices  involved 

in  this  method  may  become  quite  large.  The  largest  single  matrix  will 

be  K  ,  which  is  6M  x  6M,  and  even  for  a  system  of  only  fifteen  primary 
P 

elements,  this  matrix  would  be  180  x  180.  A  matrix  this  size  has  32,400 
elements  and  if  stored  in  the  core  memory  of  a  32K  computer  would  leave 
only  368  cells  for  other  variables,  the  monitor,  the  resident  program, 
etc.  This  is  far  too  few  cells  for  these  purposes.   If  double  precision 
arithmetic  were  to  be  used  then  twice  as  many  cells  would  be  required 
and  the  matrix  could  not  be  stored  in  core  memory. 
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CHAPTER  III 
REDUCTION  OF  SUB  SYSTEMS 

3.1  Purpose 

It  is  the  purpose  of  this  thesis  to  develop  and  exhibit  methods  by 
which  the  mathematical  manipulations  may  be  simplified  or  partitioned 
with  the  particular  purpose  of  keeping  well  within  the  core  memory  capa- 
city of  the  usual  large  engineering  computer.  Matrices  too  large  for  the 
core  can  be  handled  with  the  use  of  auxiliary  memories  or  intermediate 
tapes;  but  this  slows  computations  by  several  orders  of  magnitude  and 
is  to  be  avoided  if  possible,  particularly  in  a  problem  of  this  type 
where  the  solutions  are  gained  by  successive  iterations. 

3.2  Theory 

The  basic  propositions  which  will  be  demonstrated  are;  (1)  that  any 
ganglion  or  sub-system  can  be  represented  by  an  equivalent  stiffness  mat- 
rix of  size  6B1  x  6B1  where  B1  is  the  number  of  nodes  at  which  the  gan- 
glion is  connected  to  the  overall  system,  and  (2)  the  size  of  the  overall 
system  frequency  determinant  can  be  reduced  so  that  it  is,  at  most,  12B 
x  12B  where  B  is  the  number  of  boundary  nodes  of  the  system. 

Start  by  taking  any  ganglion  of  the  system  that  does  not  include  a 
primary  element  connected  to  a  boundary  node,  and  treat  it  as  if  it  were 
a  complete  system  having  B'  boundary  nodes,  N'  nodes,  and  M'/2  primary 
elements.   Going  through  the  manipulations  indicated  in  part  II  the 
equation  representing  this  ganglion  eventually  may  be  reduced  to  the  form: 

F  ,   -  K  .  D  .  (cf«  equation  2.7) 

N'     TJ'   H1 
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Then  by  ordering  the  matrices  so  that  the  B'  boundary  node  vectors  appear 
first  in  F  and  D  and  partitioning  we  obtain: 


NB1 


h 


11 


h 


12 


to' 


3.1.1 


*N21      KN22j   _"N2 
Now  it  must  be  understood  that  there  is  a  difference  between  this  case 
and  the  solution  of  the  overall  problem,  the  difference  being  that  while 
the  elements  of  F.TO  are  identically  equal  to  zero,  the  elements  of 
D._ ,   are  not.   In  other  words,  for  the  overall  system  the  existence  of 
a  boundary  force  indicates  a  rigid  constraint  component   (i.e.,  no  dis- 
placement) whereas  for  the  sub-system,  there  is  no  such  restriction. 
Proceeding, 


NB' 


11  DNB' 


21  DNB' 


12  DN2 


3.1.1 


DN2. 


Solving  the  second  of  these  for  D  _  we  get 

DN2  *   "*M22  *N21  DNB' 
Substituting  it  into  the  first  equation  yields 

FNB'  "  L^NU   "  ^12  *N22  h2l\       °NB' 

*  KN'   DNB'  3.1.2 

where  F   ,  and  D   ,  are  6B'  x  1  and  K^,  is  6B1  x  6B1 .   K^  is  defined  by 
this  relation. 

What  has  been  generated  might  be  thought  of  as  a  B' -ended  pseudo- 
element  to  replace  the  ganglion  and  retain  its  identical  characteris- 
tics.  The  pseudo-element  is  represented  by  a  stiffness  matrix,  but  note 
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that  this  is  no  longer  the  12  x  12  stiffness  matrix  of  a  simple  element 
but  a  6B1  x  6Bf  matrix. 

This  procedure  may  be  repeated  for  other  ganglia  which  are  separate 
from  those  considered  previously,  or  a  pseudo-element  and  its  associated 
stiffness  matrix  (or  indeed  several  such  elements)  may  be  used  as  a  build* 
ing  block  in  constructing  a  larger  ganglion.  Each  pseudo -element  calcula- 
tion eliminates  N'-B'  nodes  from  the  sub-system.  This  leads  to  the  con- 
clusion that  the  most  "efficient'1  ganglion  to  choose  is  the  one  with  the 
fewest  boundary  nodes  for  the  most  internal  nodes. 

With  the  various  ganglia  reduced  to  pseudo-elements,  it  is  then 
possible  to  proceed  as  in  Chapter  II.  A  primitive  stiffness  matrix  is 
assembled  and  pre-and  post-multiplied  by  the  topological  matrices  to 
yield  a  node  stiffness  matrix.   This  is  reordered,  transformed,  and  re- 
ordered again,  the  frequency  determinant  is  isolated,  and  the  natural 
frequencies  are  found  by  an  iterative  technique  as  described  in  section 
2.1.  What  actually  has  been  done  is  to  build  up  to  a  solution  step  by 
step,  eliminating  as  many  nodes  as  possible  along  the  way.   It  is  obvious 
that  for  the  final  step  one  can  have  the  B  boundary  nodes  of  the  system 
each  connected  by  a  primary  element  to  one  of  the  B'  boundary  nodes  (B* 4 
B)  of  a  pseudo-element  representing  the  entire  remainder  of  the  system, 
i.e.,  the  entire  system  less  the  elements  connecting  points  B  with  points 
B'. 

It  can  be  shown  that  tills  is  so  by  considering  the  system  as  a 
whole  and  then  reducing  it,  rather  than  by  starting  with  a  small  sub- 
system and  building  up  to  the  complete  system.  The  first  ganglion  can 
then  be  chosen  to  be  everything  except  those  primary  elements  that  have 
one  end  incident  on  a  boundary  node.   If  such  a  ganglion  can  not  be 
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chosen,  then  there  is  more  than  one  structurally  independent  system  in- 
volved and  they  must  be  solved  separately.   This  statement  is  not  as 
trivial  as  it  may  seem  to  be  at  first  reading.   It  is  entirely  possible 
that  a  piping  network  that  would  be  considered  as  one  system  by  reason 
of  its  application  or  fluid  flow  could  actually  be  several  independent 
structural  systems.   This  could  be  caused  by  the  pipe  being  rigidly 
anchored  at  one  or  more  points  thus  effectively  isolating  the  pipe  on 
one  side  of  the  anchor  from  the  effects  of  any  change  on  the  other  side. 

Although  many  of  the  boundary  nodes  of  a  piping  system  may  be  topolo- 
gical^ identical  (for  example  all  nodes  at  which  an  element  is  rigidly 
connected  to  "ground")  it  is  necessary  to  consider  each  as  a  separate 
entity  because  of  the  necessity  of  using  different  local  coordinate 
systems  and  then  different  transformation  matrices.  Thus  if  two  or 
more  elements  terminate  on  what  is  geometrically  the  same  boundary  node, 
they  must  be  considered  as  terminating  on  topologically  separate  nodes. 
Therefore,  there  will  be  only  one  primary  element  incident  on  each  bound- 
ary node.   (See  definition,  page  7  ).  Keeping  this  in  mind  it  can  be 
seen  that  the  number  of  nodes  connected  to  boundary  nodes  by  primary 
elements  must  be  less  than,  or  at  the  most  equal  to  B.  Thus  for  the 
reduced  system  N  is  less  than,  or  at  most  equal  to  2B.   Taking  the 
limiting  case,  N  =  2B,  the  node  stiffness  matrix  K^  will  be  12B  x  12B 
and  the  frequency  determinant  will  be  (12B-P)  x  (12B-P)  where  P  lies  in 
the  range  0  ^  P^6B.   If  all  the  boundary  nodes  represent  complete  con- 
straints in  all  degrees  of  freedom  P  will  be  6B,  the  other  limit  being 
no  constraint  on  the  boundaries  and  P  equaling  zero.   Thus  the  order  of 
the  frequency  determinant  will  be  between  6B  x  6B  and  12B  x  12B.   The 
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restriction  of  having  only  one  primary  element  incident  on  a  boundary 
node  would  not  be  necessary  for  the  solution  of  a  mathematically  similar 
problem  where  physical  location  and  metric  geometry  are  of  no  importance, 
for  example,  an  electrical  network. 

3.3  Reduction  of  Parallel  Elements 

To  provide  a  tool  for  future  computational  work,  and  to  demonstrate 
the  use  of  topological  matrices  it  will  be  shown  that  parallel  elements, 
i.e.,  primary  elements  connecting  the  same  pair  of  primary  nodes,  can  be 
reduced  to  an  equivalent  element.   The  stiffness  matrix  of  this  equi- 
valent element  will  be  the  sum  of  the  stiffness  matrices  of  the  various 
component  primary  elements.   This  is  not  a  demonstration  of  the  reduc- 
tion method  of  section  3.2  as  there  are  no  nodes  eliminated.  The  re- 
duction of  two  parallel  elements  into  one  equivalent  element  is  shown 
here.   The  logical  extension  is  then  made  that  this  equivalent  element 
can  be  taken  parallel  with  another  primary  element  and  these  may  be  re- 
duced to  a  second  equivalent  element.  As  many  parallel  elements  as 
required  can  be  handled  in  this  manner. 

Take  the  case  of  two  parallel  elements,  Fig.  3.3.1. 


Fig.  3.3.1 


1  and  2  are  nodes,  and  a,b,c,  and  d  are  element  ends. 

Z   -  U.    Zt      from  which  F  .   -  K  .  D  . 
a      ba   b  a,b      ab  a,b 

Z   -  U,   Z,      from  which  F     -  K  .  D   . 
c      dc   d  c,d      cd  c,d 

After  transforming  these  to  some  global  coordinate  system  the  primitive 

stiffness  matrix  can  be  assembled,  viz.: 
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Kabll  +  Kcdll    Kabl2  +  Kcdl2 


|_Kab21  +  Kcd21 


K  u  +  K  a 
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Since   FN  =  ^  DN 
it  has  been  shown  that 
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Thus,  the  equivalent  stiffness  matrix  for  two  parallel  elements  is  the 
sum  of  the  two  individual  stiffness  matrices.  From  the  above  we  can  see 
that  the  stiffness  matrix  for  n  parallel  elements  would  be  the  sum  of  the 
n  individual  stiffness  matrices. 

This  intuitively  attractive  result  could  probably  have  been  shown 
without  resorting  to  the  topological  matrices.  The  preceding  deriva- 
tion, however,  puts  into  a  definite  mathematical  form  what  would  have  to 
be  understood  in  another  derivation.  That  is,  the  sum  of  the  forces  on 
the  node  equals  the  sum  of  the  forces  at  the  element  ends  and  the  dis- 
placements of  all  element  ends  at  any  node  are  identically  equal. 

An  equivalent  element  with  its  own  stiffness  matrix  has  now  been 
generated,  and  may  be  used  as  one  of  the  building  blocks  of  the  system. 
However,  in  a  manner  equivalent  to  reversing  what  we  did  to  get  a  stiff- 
ness matrix  from  a  transfer  matrix  in  section  2.1  we  can  get,  if  desired, 
a  transfer  matrix  for  this  equivalent  element,  in  whichever  local  co- 
ordinate system  is  convenient.   If  the  system  is  so  configured  that  after 
the  reduction  of  the  parallel  elements  the  equivalent  element  is  then  one 
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of  two  or  more  elements  of  a  primary  path,  the  transfer  matrix  may  then 
be  used  in  the  generation  of  a  primary  element  including  the  equivalent 
element. 
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CHAPTER  IV 
COMPUTER  PROGRAM  PLANNING 

4.1   Discussion 

It  should  be  obvious  by  now  that  the  solution  of  any  but  the  simplest 
problems  by  these  methods  will  require  the  use  of  a  high  speed  digital 
computer,  and  that  the  use  of  intermediate  tapes  or  auxiliary  memory  will 
most  likely  be  required. 

The  primary  object  of  all  the  calculating  is  to  find  some  or  all  of 
the  real  values  of  frequency  which  will  make  the  frequency  determinant  go 
to  zero.   The  corresponding  "mode  shapes"  can  easily  be  computed  by  a 
repeated  pass  using  the  correct  frequency.   Since  the  computer  can  only 
calculate  with  numbers,  it  will  not  be  possible  to  set  up  an  equation  for 
the  determinant  and  solve  it  exactly.   It  will  not  even  be  possible  to  set 
up  an  equation  and  try  to  solve  it  by  substituting  very  many  values  of 
frequency.   Rather,  all  calculations  must  be  done  numerically  using  assum- 
ed values  of  frequency.   The  frequencies  which  make  the  frequency  deter- 
minant equal  to  zero  must  be  found  by  an  iterative  process.   This  means 
that  virtually  all  the  calculations  must  be  redone  for  each  assumed  fre- 
quency . 

The  writing  of  a  computer  program  to  do  the  work  involved  is  not  a 
simple  matter.   The  ideal  program  should  be  written  in  a  widely  used  en- 
gineering language  (such  as  FORTRAN)  and  arranged  to  allow  data  to  be 
presented  in  a  most  expeditious  manner.  A  good  example  of  this  type  of 
program  is  "STRESS"  by  Fenves  et.  al.  [4]  [5].   "STRESS11  solves  static 
force-displacement  problems  by  use  of  a  stiffness  matrix  method. 
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No  attempt  has  been  made  to  write  a  program  for  this  thesis.  How- 
ever, the  outline  of  such  a  program  is  shown  in  Fig.  4.1.1.  A  method  for 
systematically  approaching  this  problem  is  given  in  section  4.2. 

The  reader  interested  in  actually  programming  this  type  of  problem 
is  advised  of  the  existance  of  a  FORTRAN  program  written  by  Fink  [6]  in 
1964.   Program  VIPIPE  is  severely  restricted  in  the  size  and  topological 
configuration  of  the  systems  it  can  handle.   It  was  written  to  determine 
in-plane  and  out-of -plane  vibrations  of  a  two  dimensional  piping  system 
and  does  not  make  use  of  topological  matrices.   It  would  be  a  good  start- 
ing point  for  the  generation  of  a  more  sophisticated  program  and  many  of 
its  subroutines,  in  particular  those  for  calculating  transfer  matrices, 
could  probably  be  put  to  good  use. 
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Figure  4.1.1 
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4. 2  Expansion  of  Computational  Concepts 

The  purpose  of  this  section  is  to  develop  to  a  greater  extent  the 
approach  to  a  computer  solution  for  a  general  problem  of  this  category. 
The  basic  flow  diagram  for  a  computer  program  shown  in  Fig.  4.1.1  will 
be  used  as  a  rough  guide. 

The  following  assumptions  are  made: 

a.  A  computer  memory  of  about  30K  available  cells  is  to  be 
utilized. 

b.  Intermediate  tapes  rather  than  drums  or  discs  will  be 
used.   This  assumption  puts  some  constraints  on  the 
intermediate  storage  of  information  as  binary  tapes  cannot 
be  randomly  accessed  whereas  drums  or  discs  can. 

c.  Most  of  the  calculating  will  be  done  with  double  precision 
arithmetic. 

d.  The  system  to  be  solved  may  be  of  any  size. 
The  following  systematic  approach  is  recommended. 

1.)  Display  the  system  topologically  showing  only  the  boundary 
and  primary  nodes  and  the  primary  elements  connecting  them.   Number  the 
nodes,  using  the  lowest  numbers  for  the  boundary  nodes,  and  number  the 
primary  elements. 

2.)  Determine  the  way  in  which  the  system  will  be  divided  into 
ganglia  and  subsequently  recombined.  A  procedure  for  making  cuts  for 
an  "optimum"  reduction  is  beyond  the  scope  of  this  work,  and  indeed  per- 
haps the  idea  of  an  optimal  reduction  is  frivolous.   However,  certain 
operations  should  be  readily  apparent,  and  the  first  of  these  is  the 
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combining  of  parallel  elements,  where  present.  From  this  point  forward 
it  is  an  empirical  process.  One  must  use  his  judgement  in  working  to- 
ward replacing  the  entire  system,  except  for  the  boundary  nodes  and  the 
primary  elements  incident  on  the  boundary  nodes,  with  a  single  pseudo- 
element.  For  reasons  that  will  be  discussed  later,  it  is  best  to  set 
up  two  or  more  different  sets  of  ganglia  so  that  the  problem  may  be  solved 
more  than  once  and  the  answers  compared. 

3.)  Display  the  system  so  that  the  geometry  will  be  readily  avail- 
able for  computer  programming.  An  isometric  line  drawing,  (not  neces- 
sarily to  scale,  but  at  least  reasonably  proportional  to  the  actual  layout) 
will  probably  be  adequate. 

4.)  Analyze  the  boundary  node  constraints.   It  is  necessary  that 
each  node  be  either  unrestrained  or  rigidly  restrained  in  each  of  its 
six  degrees  of  freedom.   If  this  is  not  the  case  then  an  imaginary  ele- 
ment with  characteristics  equivalent  to  the  desired  boundary  constraints 
must  be  added  between  the  element  end  and  the  boundary  node  to  satisfy 
the  above  condition. 

5.)  Analyze  each  other  node.  The  assumption  made  in  this  method 
is  that  elements  are  rigidly  connected  together  at  the  nodes.   If  this 
is  not  the  case,  an  imaginary  spring  element  of  desired  characteristics 
must  be  included  wh«r«rer  required.* 

6.)  Establish  an  origin  of  coordinates  and  an  orthogonal  set  of 
global  coordinates. 


*This  is  related  to  the  difficult  problem  of  releases.  See  reference 
[7]. 
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7.)  Determine  and  tabulate  for  each  element  the  following  data: 

a.  Type  of  element 

1.  straight  pipe 

2.  arc  of  circle  pipe 

3.  "imaginary"  element 

4.  non-pipe  (hanger) 

b.  weight  per  unit  length  of  pipe 

c.  weight  per  unit  length  of  fluid  and  insulation 

d.  pipe  diameter 

e.  wall  thickness 

f.  E  (modulus  of  elasticity) 

g.  G  (modulus  of  elasticity  in  shear) 

h.   coordinates  of  ends  and  of  center  of  curvature  (if  an 

arc) 
i.   radius  of  curvature 

8.)  Read  in  the  data  for  each  element;  make  whatever  calculations 
are  necessary  to  convert  the  raw  data  into  an  efficient  form  for  future 
calculating;  and  store  all  this  information  on  tape.  All  this  data  is 
independent  of  frequency,  and  must  be  used  in  each  iteration  to  calcu- 
late transfer  matrices.   It  is  thus  vital  that  the  information  be  stored 
in  the  order  that  it  is  to  be  used  so  that  the  tape  may  be  used  in  the 
most  expeditious  manner.   Let  us  assume  that  we  have  a  set  of  data  cards 
carrying  all  the  information  for  all  the  elements  of  each  primary  ele- 
ment. We  also  have  a  larger  set  composed  of  all  these  sets,  so  that  the 
data  is  on  cards  for  all  elements  of  the  system.   Now  go  back  to  2)  and 
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pick  the  first  ganglion  to  be  reduced.  Decide  on  the  order  in  which  the 
primitive  stiffness  matrix  is  to  be  assembled.   (This  will  also  be  a 
good  time  to  determine  the  way  in  which  the  node  stiffness  matrix  is  to 
be  ordered  and  to  generate  the  topological  matrices  for  this  ganglion, 
though  these  items  will  not  be  needed  for  a  while.)  Assemble  the  data 
cards  for  each  primary  element  in  order  from  one  end  to  the  other. 
Remember  that  the  coordinates  of  the  first  end  of  the  first  element  will 
determine  the  local  coordinate  system  in  which  that  particular  primary 
element  will  be  calculated.  For  each  primary  element  set  there  should 
be  an  initial  card  indicating  the  numbering   jf  elements  in  the  primary 
element.  Assemble  all  the  primary  element  card  sets  in  the  order  in 
which  the  primitive  stiffness  matrix  is  to  be  assembled. 

9.)  Choose  the  next  reduction  to  be  made.  If  it  is  independent  of 
the  primary  elements  already  assembled  for  reading,  proceed  as  above.  If 
it  contains  the  ganglion  (or  several  ganglia)  assembled  for  reading, 
proceed  as  above  only  for  the  primary  elements  not  already  used.   If  it 
is  composed  entirely  of  elements  whose  cards  are  already  assembled,  no 
action  is  required  at  this  time.   (  It  would  still  be  good  at  this  time 
to  set  up  the  order  for  the  node  stiffness  matrix  and  generate  the  topo- 
logical matrices  for  this  ganglion) . 

10.)  When  the  ordering  is  complete  the  element  data  is  ready  to  be 
read  in.  A  flag  indicating  that  no  more  element  data  is  to  be  read  should 
be  added.  The  number  of  elements  in  each  primary  element  are  stored  in 
an  array  to  facilitate  retrieval  of  information  from  the  tape.  The  decis- 
ion on  exactly  how  to  handle  data,  particularly  in  reference  to  the  amount 
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of  data  to  be  read  in  and  transferred  to  tape  at  one  time  must  be  left  to 
the  programmer.   It  is  faster  to  accumulate  data  in  core  and  transfer  it 
all  to  tape  with  one  WRITE  TAPE  statement  than  to  go  to  tape  with  each 
individual  set  as  it  is  calculated.   It  should  be  noted  that  much  economy 
in  data  card  preparation  can  be  obtained  by  eliminating  the  necessity  for 
duplicating  information  that  is  identical  for  many  elements,  i.e.,  Young's 
modulus,  diameter,  etc. 

11.)  Now  starting  with  the  first  ganglion  to  be  reduced,  read  in 

T 
the  appropriate  topological  matrix  (A  or  A  )  and  store  it.   It  will  not 

T 
be  necessary  to  have  both  A  and  A  in  memory  since  by  programming  we  can 

utilize  the  stored  matrix  either  as  itself  or  as  its  transpose.   Since  all 

elements  of  the  topological  matrices  are  either  61  or  60   great  storage 

economy  can  be  obtained  by  representing  these  6x6  elements  by  a  Boolean 

1  or  0  respectively  and  storing  thirty- two  of  these  per  core  memory  word.* 

It  will  probably  not  be  necessary  to  go  out  on  tape  with  these  matrices. 

At  this  point  all  the  data  that  will  be  constant  from  iteration  to 

iteration  should  be  either  on  tape  or  in  core. 

12.)  With  the  assumption  of  an  initial  frequency  the  first  iteration 
may  now  be  started.   From  this  point  on,  all  calculations  will  be  done 
with  double  precision  arithmetic.   Start  by  taking  the  first  ganglion  as 
previously  ordered,  and  read  back  into  the  core  all  the  data  relevant  to 
the  elements  of  this  ganglion.  By  the  use  of  sub-routines  the  stiffness 
matrix  for  each  primary  element  of  the  ganglion  is  calculated.   It  will 


*The  number  of  logical  elements  that  may  be  stored  in  one  word  varies 
from  computer  to  computer.  The  number  thirty- two  used  here  applies  to 
the  CDC  1604. 
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not  be  necessary  to  store  the  numerous  transfer,  rotation  and  transforma- 
tion matrices  that  are  generated  for  if  they  are  generated  in  the  correct 
order  and  multiplied  together  then  all  that  need  be  stored  is  the  pro- 
duct. The  stiffness  matrix  for  each  primary  element  will  occupy  2  x  12 
x  12  =  288  core  storage  locations  where  the  2  results  from  double  pre- 
cision arithmetic.  Mow  although  the  primitive  stiffness  matrix,  K  ,  is 
shown  as  a  large  square  matrix,  it  will  be  assembled  in  the  computer  as 
a  three  dimensional  array.  The  reason  for  this  is  that  K  consists  only 
of  diagonal  elements,  each  of  which  is  of  order  12  x  12,  and  all  the  other 
elements  are  zero.  By  using  a  three  dimensional  array  we  can  store  only 
the  diagonal  elements  and  not  waste  space  storing  zeros.  The  storage  re- 
quired for  K  will  be  288  x  M'/2,  and  if  the  stiffness  matrices  have 
been  cleverly  stored  when  they  were  generated  K  already  exists.  The 
node  stiffness  matrix  is  now  generated  from  K^  -     A  K  A.   Note  that  the 
column  vectors  of  node  forces  and  displacements  are  never  actually  used 


in  the  solution.  From  K^  we  generate  K^,   as  shown  in  equations  3.1.1 
and  3.1.2.  The  order  of  this  matrix  will  be  6B'  x  6B';  and  the  number  of 
memory  words  necessary  to  hold  it  will  be  2  x  6B'  x  6B'.  As  before,  the 
2  represents  the  requirement  of  two  words  per  item  in  double  precision 
arithmetic. 

Assuming  that  this  reduced  ganglion  will  be  used  as  part  of  the 
next  ganglion  to  be  reduced,  we  will  plan  to  leave  K..,   in  core  memory. 
If  it  is  not  to  be  used  for  a  while,  and  if  space  is  at  a  premium  for 
other  calculations,  it  may  be  better  to  store  K  ,   on  tape. 

13.)  By  continuing  in  a  similar  manner  eventually  we  will  arrive 
at  the  largest  ganglion,  namely,  that  one  including  everything  except 
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the  primary  elements  incident  upon  the  boundary  nodes.   The  order  of  the 
resultant  stiffness  matrix  for  this  ganglion  will  be  6B  x  6B.   Using 
this  stiffness  matrix,  and  the  stiffness  matrices  of  the  remaining 
primary  elements  we  now  set  up  the  primitive  stiffness  matrix  for  the 
entire  system  and  by  pre-  and  post-multiplying  it  by  the  topological 
matrices  obtain  the  node  stiffness  matrix  for  the  entire  system.   Since 
we  have  assembled  these  matrices  carefully,  the  reordering  described  in 
section  2.h  will  not  be  necessary.  The  transformation  back  to  local  co- 
ordinates is  made  for  the  boundary  nodes.  With  the  boundary  nodes  re- 
presented in  local  coordinates  it  is  now  possible  to  make  the  second  re- 
ordering as  described  in  section  2.i.  The  frequency  determinant  is  now 
isolated  and  its  value  calculated.  This  value  and  the  corresponding 
value  of  frequency  at  which  the  calculations  have  been  made  are  saved. 

14.)  A  new  value  of  frequency  is  now  assumed  and  the  above  process 
repeated.   Subsequent  values  of  the  determinant  may  be  used  in  conjunc- 
tion with  the  values  found  previously  to  predict  the  next  frequency  to  be 
tried  in  order  that  the  calculations  may  converge  on  the  desired  natural 
frequency.  Alternatively,  calculations  may  be  made  at  constant  intervals 
and  the  results  plotted. 

4.3  Computational  Difficulties 

There  is  one  predominant  problem  that  keeps  appearing,  if  only 
implicitly,  in  work  with  matrices.  Any  time  matrix  equations  are  mani- 
pulated in  such  a  way  that  the  inverse  of  some  matrix  is  required,  the 
assumption  must  be  made  that  the  matrix  to  be  inverted  is  non-singular. 
Matrices  are  inverted  many  times  in  the  course  of  a  solution  by  the  methods 
discussed  herein,  and  the  possibilities  of  having  singular  matrices  occur 
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occasionally  are  significant  and  somewhat  disturbing. 

In  the  determination  of  the  stiffness  matrix  of  an  element  from  its 
transfer  matrix  it  is  necessary  to  invert  a  6  x  6  matrix  as  described  in 
section  2.b. 
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Note  that  \J.„   is  the  6x6  matrix  to  be  inverted.  As  some  or  all  of  the 
elements  of  IL  _  are  frequency  dependent  (all  other  parameters  are  depen- 
dent upon  physical  constants  and  are  thus  fixed)  it  can  be  seen  that  there 
will  be  some  values  of  frequency  for  which  the  determinant  of  U.„  will 
vanish.   These  frequencies  correspond  to  the  natural  frequencies  of  the 
element  with  both  ends  "built-in"  or  cantilevered.  In  matrix  form,  from 
2.2- 
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■  0.  Only  the  real  values  of  frequency  that  satisfy 
this  relationship  are  of  interest. 

What  this  means  in  terms  of  programming  is  that  whenever  the  assum- 
ed value  of  frequency  at  which  a  series  of  calculations  is  being  made  is 
the  same  as  (or  close  to)  one  of  the  natural  frequencies  of  one  of  the 
primary  elements  it  will  not  be  possible  to  generate  the  stiffness  matrix 
for  that  element.   It  will  be  necessary  to  assume  a  new  frequency  near  the 
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unacceptable  one  and  to  try  to  calculate  to  a  conclusion  again. 

A  similar  problem  may  arise  in  the  calculation  of  the  equivalent 
stiffness  matrix  for  a  ganglion.  Assuming  that  stiffness  matrices  have 
been  satisfactory  calculated  for  each  primary  element  of  the  ganglion 
as  shown  in  the  arithmetical  manipulations  of  equations  3.1.1,  it  is 
still  necessary  to  invert  a  matrix  and  this  matrix,  also  being  made  up  of 
elements  dependent  upon  the  frequency,  may  become  singular.   The  fre- 
quencies at  which  this  matrix  becomes  singular  correspond  to  the  natural 
frequencies  of  a  system  equivalent  to  the  ganglion  but  with  the  boundary 
nodes  rigidly  fixed. 

The  inability  to  calculate  to  a  conclusion  due  to  an  inability  to 
invert  part  of  the  node  stiffness  matrix  of  a  ganglion  can  be  avoided, 
in  most  cases,  by  breaking  the  system  up  differently  and  thus  using  dif- 
ferent ganglia.   It  is  for  this  reason  that  in  section  4.2  it  is  recommend- 
ed that  any  problem  of  this  type  be  solved  using  several  different  combina- 
tions of  ganglia.   Any  reasonably  sophisticated  computer  program  will,  of 
course,  provide  some  sort  of  diagnostic  message  and  proceed  to  calculate 
at  a  different  value  of  frequency  whenever  some  matrix  cannot  be  inverted. 

The  problem  of  a  singular  matrix  affecting  the  calculations  should 
not  be  too  serious  at  the  low  frequency  end  of  the  vibrational  spectrum. 
It  can  reasonably  be  expected  that  the  natural  frequencies  of  a  ganglion 
will  be  considerably  greater  than  the  lowest  natural  frequency  of  the 
overall  system.   It  is  not  possible  to  forecast  at  what  frequency  pro- 
blems will  actually  begin  to  occur.  What  is  to  be  done  is  just  to  pro- 
ceed with  the  calculations,  avoiding  those  frequencies  at  which  calcula- 
tions cannot  be  made,  but  calculating  as  closely  as  possible  on  both 
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sides  of  such  frequencies. 

There  are  other  problems  involved  in  the  creation  of  a  computer  pro- 
gram capable  of  handling  large  linear  systems  by  matrix  methods.  There 
is,  for  example,  the  possibility  that  even  after  all  reductions  have  been 
made  the  frequency  determinant  will  still  be  too  large  for  core  memory. 
There  will  also  be  the  need  for  making  a  number  of  tests  during  the  course 
of  a  program  run,  both  for  debugging  purposes  and  to  provide  a  hedge 
against  gross  errors.   (An  example  of  this  would  be  the  testing  of  the 
value  of  mass  of  an  element  to  insure  that  it  is  positive).   It  is  re- 
commended therefore  that  the  "team"  for  developing  this  program  consist 
of  at  least  a  mechanical  engineer  well  versed  in  vibrations  and  a  com- 
puter programmer  with  experience  in  handling  large  systems  and  matrices. 
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